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The completion of a network of advanced laser-interferometric gravitational-wave obser- 
vatories around 2001 will make possible the study of the inspiral and coalescence of binary 
^ ' systems of compact objects (neutron stars and black holes), using gravitational radiation. 

To extract useful information from the waves, such as the masses and spins of the bodies, 
theoretical general relativistic gravitational waveform templates of extremely high accuracy 
1/^ ' will be needed for filtering the data, probably as accurate as 0[(«/c)^] beyond the predictions 

of the quadrupole formula. We summarize a method, called DIRE, for Direct Integration 
of the Relaxed Einstein Equations, which extends and improves an earlier framework due 
to Epstein and Wagoner, in which Einstein's equations are recast as a flat spacetime wave 
^ ' equation with source composed of matter confined to compact regions and gravitational 

C""~~- \ non-linearities extending to infinity. The new method is free of divergences or undefined 

VO . integrals, correctly predicts all gravitational wave "tail" efi^ects caused by backscatter of the 

^3 ' outgoing radiation off the background curved spacetime, and yields radiation that propa- 

^^ ' gates asymptotically along true null cones of the curved spacetime. The method also yields 

equations of motion through 0[{v/c)'^\, radiation-reaction terms at 0[{v/c)^] and 0[{v/cf\, 
0^ ■ and gravitational waveforms and energy fiux through 0[(ii/c)*], in agreement with other 

CTN ' approaches. We report on progress in evaluating the 0[{v/c)^] contributions. 

o 



§1. Introduction 



Some time in the next decade, a new window for astronomy and relativistic 
. gravity may be realized, with the completion and operation of kilometer-scale, laser 

r> \ interferometric gravitational-wave observatories in the U.S. (LIGO project), Europe 

C3 : (VIRGO and GEO-600 projects) and Japan (TAMA 300 project). Gravitational- 

wave searches at these observatories are scheduled to commence around 2002. The 
LIGO broad-band antennae will have the capability of detecting and measuring 
the gravitational waveforms from astronomical sources in a frequency band between 
about 10 Hz (the seismic noise cutoff) and 500 Hz (the photon counting noise cutoff), 
with a maximum sensitivity to strain at around 100 Hz of Al/l ~ 10~^^ (rms). The 
most promising source for detection and study of the gravitational- wave signal is the 
"inspiralling compact binary" - a binary system of neutron stars or black holes (or 
one of each) in the final minutes of a death dance leading to a violent merger. Such 
is the fate, for example of the Hulse- Taylor binary pulsar PSR 1913-1-16 in about 240 
million years. Given the expected sensitivity of the "advanced LIGO" (around 2007), 
which could see such sources out to hundreds of megaparsecs, it has been estimated 
that from 3 to 100 annual inspiral events could be detectable. Other sources, such 
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as supernova core collapse events, instabilities in rapidly rotating nascent neutron 
stars, signals from non-axisymmetric pulsars, and a stochastic background of waves, 
may be detectable (for reviews, see Ref. |I|) and other articles in this volume). 

The analysis of gravitational-wave data from such inspiral sources will involve 
some form of matched filtering of the noisy detector output against an ensemble of 
theoretical "template" waveforms which depend on the intrinsic parameters of the 
inspiralling binary, such as the component masses, spins, and so on, and on its inspiral 
evolution. How accurate must a template be in order to "match" the waveform from 
a given source (where by a match we mean maximizing the cross-correlation or the 
signal-to-noise ratio)? In the total accumulated phase of the wave detected in the 
sensitive bandwidth, the template must match the signal to a fraction of a cycle. For 
two inspiralling neutron stars, around 16,000 cycles should be detected; this implies 
a phasing accuracy of 10~^ or better. Since v/c ~ 1/10 during the late inspiral, 
this means that correction terms in the phasing at the level of (v/c)^ or higher are 
needed. More formal analyses confirm this intuition. B^ 

Because it is a slow-motion system {v/c ~ 10~^), the binary pulsar is sensi- 
tive only to the lowest-order effects of gravitational radiation as predicted by the 
quadrupole formula. Nevertheless, the first correction terms £if order v/c and {v/c)"^ 
to the quadrupole formula, were calculated as early as 19760^. These are now con- 
ventionally called "post-Newtonian" (PN) corrections, with each power of v/c corre- 
sponding to half a post-Newtonian order (0.5PN), in analogy with post-Newtonian 
corrections to the Newtonian equations of motion. eP In 1976, the post-Newtonian 
corrections to the quadrupole formula were of purely academic, rather than obser- 
vational interest. 

But for laser-interferometric observations of gravitational waves, the bottom line 
is that, in order to measure the astrophysical parameters of the source and to test the 
properties of the gravitational waves, it is necessary to derive the gravitational wave- 
form and the resulting radiation back-reaction on the orbit phasing at least to 2PN, 
or second post-Newtonian order, 0[(i'/c)^], beyond the quadrupole approximation, 
and probably to 3PN order. 

§2. Post-Newtonian Generation of Gravitational Waves 

The motion of isolated binary systems and the generation of gravitational ra- 
diation are long-standing problems that date back to the first years following the 
publication of GR, when Einstein calculated the gravitational radiation emitted by 
a laboratory-scale object using the linearized version of GR. Shortly after the dis- 
covery of the binary pulsar PSR 1913+16 in 1974, questions were raised about the 
foundations of the "quadrupole formula" for gravitational radiation damping (and in 
some quarters, even about its quantitative validity). These questions were answered 
in part by theoretical work designed to shore up the foundations of the quadrupole 
approximation, EP and in part (perhaps mostly) by the agreement between the pre- 
dictions of the quadrupole formula and the observed rate of damping of the pulsar's 
orbit. 

The challenge of providing accurate templates for LIGO- VIRGO data analysis 
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has led to major efforts to calculate gravitational waves to high PN order. Three 
approaches have been developed. 

The BDI approach of Blanchet, Damour and Iyer is based on a mixed post- 
Newtonian and "post-Minkowskian" framework for solving Einstein's equations ap- 
proximately, developed in a series of papers by Damour and colleagues. Q) The idea 
is to solve the vacuum Einstein equations in the exterior of the material sources 
extending out to the radiation zone in an expansion ("post-Minkowskian") in "non- 
linearity" (effectively an expansion in powers of Newton's constant G), and to express 
the asymptotic solutions in terms of a set of formal, time-dependent, symmetric and 
trace- free (STF) multipole moments. LP Then, in a near zone within one character- 
istic wavelength of the radiation, the equations including the material source are 
solved in a slow-motion approximation (expansion in powers of 1/c) that yields a set 
of STF source multipole moments expressed as integrals over the "effective" source, 
including both matter and gravitational field contributions. The solutions involving 
the two sets of moments are then matched in an intermediate zone, resulting in a 
connection between the formal radiative moments and the source moments. The 
matching also provides a natural way, using analytic continuation, to regularize in- 
tegrals involving the non-compact contributions of gravitational stress-energy, that 
might otherwise be divergent. For further details, see the article by Blanchet in this 
volume. 

An approach called DIRE is based on a framework developed by Epstein and 
Wagoner (EW)B\ and extended by Will, Wiseman and Pati. We shall describe DIRE 
briefly below. 

A third approach, valid only in the limit in which one mass is much smaller than 
the other, is that of black-hole perturbation theory. This method provides numerical 
results that are exact in v/c, as well as analytical results expressed as series in powers 
of v/c, both for non-rotating and for rotating black holes. For non-rotating holes, 
the analytical expansions have been carried to 5.5PN order. EP In all cases of suitable 
overlap, the results of all three methods agree precisely. 

§3. Direct Integration of the Relaxed Einstein Equations (DIRE) 

Like the BDI approach, DIRE involves rewriting the Einstein equations in their 
"relaxed" form, namely as an inhomogeneous, flat-spacetime wave equation for a 
field /i"^, whose formal solution can be written 

raP(+ - Iy _ y'I yM 



h'^^it, x) =4 / ^° (t 1^ 7^1 '^'^ ^3^^ , (3-1) 



c 

where the source r"^ consists of both the material stress-energy, and a "gravitational 
stress-energy" made up of all the terms non- linear in /i"^, and the integration is over 
the past flat-spacetime null cone C of the field point (t,x) (see Fig. 1). The wave 
equation is accompanied by a harmonic or deDonder gauge condition /i"^,/3 = 0, 
which serves to specify a coordinate system, and also imposes equations of motion 
on the sources. Unlike the BDI approach, a single formal solution is written down, 
valid everywhere in spacetime. This formal solution, is then iterated in a slow- 
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motion {v/c < 1), weak-field (||/i"^|| < 1 ) approximation, that is very similar to the 
corresponding procedure in electromagnetism. However, because the integrand of 
this retarded integral is not compact by virtue of the non-linear field contributions, 
the original EW formalism quickly runs up against integrals that are not well defined, 
or worse, are divergent. Although at the lowest quadrupole and first few PN orders, 
various arguments can be given to justify sweeping such problems under the rug, eP 
they are not very rigorous, and provide no guarantee that the divergences do not 
become insurmountable at higher orders. As a consequence, despite efforts to cure 
the problem, the EW formalism fell into some disfavor as a route to higher orders, 
although an extension to 1.5PN order was accomplished. 113* 

The resolution of this problem involves taking literally the statement that the 
solution is Kjretarded integral, i.e. an integral over the entire past null cone of the 
field point. EiP To be sure, that part of the integral that extends over the intersection 
M between the past null cone and the material source and the near zone is still 
approximated as usual by a slow-motion expansion involving spatial integrals over a 
constant-time hypersurface M of moments of the source, including the non-compact 
gravitational contributions, just as in the BDI framework (Fig. 2). But instead of 
cavalierly extending the spatial integrals to infinity as was implicit in the original EW 
framework, and risking undefined or divergent integrals, we terminate the integrals 
at the boundary of the near zone, chosen to be at a radius TZ given roughly by 
one wavelength of the gravitational radiation. For the integral over the rest of the 
past null cone exterior to the near zone ("radiation zone"), we neither make a slow- 
motion expansion nor continue to integrate over a spatial hypersurface, instead we 
use a coordinate transformation in the integral from the spatial coordinates d^x' to 
quasi-null coordinates du' , d6', dcj)' , where 

t-n' = r' + |x-x'|, (3-2) 

to convert the integral into a convenient, easy-to-calculate form, that is manifestly 
convergent, subject only to reasonable assumptions about the past behavior of the 
source: 

hf_j^{t, x)=4 r du' I ^"^(^^' + ^^^0 [/(^'^ Q>)fd^Q' . (3-3) 

This transformation was suggested by earlier work on a non-linear gravitational-wave 
phenomenon called the Christodoulou memory. Y3 Not only are all integrations now 
explicitly finite and convergent, one can show that all contributions from the finite, 
near-zone spatial integrals that depend upon TZ are actually cancelled by correspond- 
ing terms from the radiation-zone integrals, valid for both positive and negative 
powers of TZ and for terms logarithmic in TZ. t3^ Thus the procedure, as expected, 
has no dependence on the artificially chosen boundary radius TZ of the near-zone. In 
addition, the method can be carried to higher orders in a straightforward manner. 
The result is a manifestly finite, well-defined procedure for calculating gravitational 
radiation to high orders. 

Thus, for field points in the far zone, the integral over the near zone takes the 
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Field point (t,x) 




Matter source 

Fig. 1. Past harmonic null cone C of the field point (t,x) intersects the near zone T) in the hyper- 
surface jV. 

(t,x) 





^ 



t-r 



Fig. 2. Taylor expansion of retarded time dependence on M results in multipole moments integrated 
over the spatial hypersurface M, corresponding to fixed retarded time t — r 
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standard form of a multipole expansion, 



oo / 1 "\ m am / -i \ 

"^^C.-) = 4E^ 3,...„a,.. (;«"».-".(, - d) , (3.4a) 

'(t-r)=/ r"^(t-r,x')2;'=i...a;'=™dV, (3-4b) 



m=0 

jiral3ki...k„ 



where integrals are over the finite hypersurface A4. For field points in the near 
zone, the integral over the near zone can be expanded in the form of a sequence of 
"Poisson"-like potentials and superpotentials (and "superduper" -potentials) evalu- 
ated at a fixed time t: 



^^'(*'-) = 4E ^^/^-"'(^'-Olx-xr-^dV. (3-5) 

m=0 "^' ^* ^-^ 

In each case, the integrals must be combined with the corresponding integral over 
the rest of the past light cone, Eq. (^]^). Because of the aforementioned general 
proof, it is not necessary to keep any terms in these integrals that depend explicitly 
on the radius 7^; this simplifies calculations considerably. 

§4. Equations of Motion to 3.5PN Order 

We assume that the orbiting bodies are sufficiently small compared to their 
separation that tidal effects, or effects due to their finite size, can be ignored. For 
inspiralling compact binaries, this is believed to be a good approximation until the 
final few orbits. This amounts to replacing a perfect fluid stress-energy tensor, 

T^"^ = {p + p)u^u'' +pg^\ (4-1) 

with that of "point-masses" : 

r^'^ = ^mA<5^(x-XA)'u^n7'u°^/^. (4-2) 

A 

However, because of gravitational non-linearities, such a stress-energy tensor will lead 
to infinities at the location of each body, hence one must find a way to regularize in 
order to isolate the physically relevant terms. Blanchet et al. use a regularization 
procedure based on the Hadamard "partie fini" . Our approach, which is less formal, 
though probably equivalent, is to isolate those terms in any integral of fields over 
a body that neither vanish nor blow up as D, the size of the body, shrinks to 
zero. Terms that vanish as D represent tidal and spin effects (and their relativistic 
generalizations), which we are ignoring. Terms that diverge as D~ are "self-energy" 
terms; we assume that these can be uniformly absorbed into renormalized masses for 
the bodies. It is important to stress that this is an assumption, whose validity has 
been checked in general only to IPN order (no Nordtvedt effect in OR) and under 
restricted circumstances to 2PN order. The result is a well-defined procedure for 
keeping "finite, point-mass" terms. 
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With these assumptions, the equations of motion for each body take the form 
of a geodesic equation, 






+ p*ri,v'^v^-p*rlv^v'^v'' = o 







(4-3) 



where p* = pil^^/^ -^ YliA "^Al5^(x — x^). 

To obtain equations of motion vahd through 3.5PN order, it is necessary to 
iterate the relaxed Einstein equation four times. Evaluating the resulting Poisson- 



like potentials for two fluid balls, integrating the equation of motion (4-3) over one 
of the bodies, and keeping only terms that are finite as the bodies shrink in size, one 
obtains equations of motion of the schematic form 



771,2 



^[n* + 0(e) + O(e^) + Oii'l') + O(e^) + 0[e'l') + . . .] ,(4-4a) 

(4-4b) 



l^[l-21 



Uc 



m/r 



x'2)/r. The expansion 



V is the relative 



where r is the distance between the bodies and n* 
parameter e is related to the orbital variables by e 
velocity, and m = mi + m2 is the total mass (G = c = 1). 

We have evaluated all contributions to h^'^ formally through 3.5PN order in 
terms of Poisson-like potentials, and have calculated them explicitly for two compact 
bodies through 2.5PN order and at 3.5PN order. The formidable task of evaluating 
the contributions at 3PN order is in progress. At 2PN order, we obtain equations 
of motion in complete agreement with those of Damour and Deruelle (Eqs. (154) - 
(160) of Ref. Il)) and Blanchet et al. (Eq. (8.4) of Ref. H)). 

The contributions at 2.5PN and 3.5PN order representeravitional-radiation re- 
action and its post-Newtonian corrections. Iyer and WillH^ have shown that, as- 
suming energy and angular momentum balance, the relative two-body equations of 
motion at 2.5PN order can be written in the form 



8 
a = --r]{m/r'^){m/r) -{A2.5 + A3,5)rn + {B2.5 + -B3.5)v 

5 L 



where rj = mi'm2/m^, r = dr/dt, and 



A. 



2.5 



23 m 

(3 + Spy + (— + 2a - 3/?) 5Pr^ , 

3 r 



m 



B2.5 = {2 + aX + {2-a) 3(l + Q)?^^ 

r 



(4.5) 

(4- 6a) 
(4-6b) 



where a and /3 are arbitrary, and reflect the effects of coordinate freedom on the equa- 
tions of motion. The values a = 4, /3 = 5 correspond to the so-called "Burke-Thorne" 
gauge, in which the radiation reaction is expressed solely as a quasi-Newtonian po- 
tential 'pRR = — gM^- (t)x*x-^, where Mij{t) is the traceless moment of inertia tensor 
of the system and the superscript (5) denotes five time derivatives. This also corre- 
sponds to the gauge used by Blanchet O-' . Our 2.5PN equations of motion yield the 
values a = —1, /3 = 0, which corresponds to the gauge used by Damour and Deruelle 
(Eq. (161) of Ref. |l|)). 
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At 3.5PN order, the expressions for A3. 5 and -B3.5 are 

^3.5 = aiv'^ + a2V^rn/r + a^v^r"^ + a^r'^Tn/r + a^^r'^ + aQ{'m/r)'^ ,{4:-7a) 
B3.5 = biv^ + b2v'^'m/r + b3v'^r'^ + b4r'^'m/r + b^r'^ + be{m/r)^ . (4-7b) 

Energy and angular momentum balance yield values for the 12 coefficients modulo 
6 arbitrary gauge parameters. Our 3.5PN equations of motion yield the values 

183 15 , 313 3 ,, , 

"^ = -^-y^' ^^ = -^-2^' ('-'"^ 

173 186 , 205 37 , , , , 

a2 = -^-—V, b2 = — + -v, (4-8b) 

285 15 , 339 3 ,, , 

as = — + y^' ^3 = — +2^' (4-8c) 

147 _ , 205 106 ,, ,, 

«4 = -—-47??, ^4 = -^-—??, (4-8d) 

as = -70, 65 = -75, (4-8e) 

«6 = -^-23r?, 66 = -^-13r?. (4-81) 

By comparing these values with the Iyer and Will expressions (Eqs. (2.18) of Ref. 
[l6|)), we obtain values for the 6 arbitrary gauge parameters. The fact that 12 con- 
straints yield a consistent solution for the 6 parameters is a useful check of the method 
and algebra. The result yields radiation reaction equations to 3.5PN or der in the 
generalization of Damour-Deruelle gauge. Blanchet's multipole expressions IIZP for ra- 
diation reaction at 3.5PN order yield a different set of coefficients, which correspond 
to the generalization of Burke-Thorne gauge. 

§5. Gravitational Radiation Waveform and Energy Flux 

By iterating the relaxed Einstein equations for field points in the far zone, and 
substituting the two-body equations of motion to the appropriate order, we obtain 
an explicit formula for the transverse-traceless (TT) part of the radiation-zone field, 
denoted /i*-', which is the waveform to be detected in laser interferometric systems. 
In terms of an expansion beyond the quadrupole formula, it has the schematic form, 

h'^ = ^{Q'^[l + Oie'/')+Oie) + Oie'/') + Oie')...]}^^, (5-1) 

where fi is the reduced mass, R is the distance to the source, and Q^^ represents 
two time derivatives of the mass quadrupole moment tensor (the series actually 
contains multipole orders beyond quadrupole). The 0.5PN and IPN terms were de- 
rived by Wagoner and Will,0' the 1.5PN terms by Wiseman. t3) The contribution of 
gravitational- wave "tails" , caused by backscatter of the outgoing radiation off the 
background spacetime curvature, at 0(e'^'^), were derived and studied by several 
authors. The 2PN terms including 2PN tail contributions were derived by two in- 



dependent groups and are in complete agreement (see Refs. 11) and O) for details 
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and references to earlier work). The 2.5PN terms and various specific higher-order 
terms, such as "tails-of-tails" have been derived by Blanchet and collaborators, while 
the formidable job of evaluating the 3PN terms is still in progress. 

To illustrate some of the results of the DIRE method, we note that, in calculating 
the gravitational waveform through 2PN order, we explicitly retained all terms that 
depend on positive powers of the radius 7^ of the near zone.li^ The multipole 
expansion of the integral over the near zone yielded 



/i^(t, x) = 7^— independent terms + 1/7^ terms 



1912 m 
315 'R 



(4)n*i 



Q'^{u)n, 



(5-2) 



where ^'^'Q^^{u) represents four time derivative of the quadrupole moment of the 
source. The 7?.-dependent term is of 2PN order; this explains why sweeping "diver- 
gent" terms under the rug at IPN order was successful, if not fully justified. The 
integral over the rest of the null cone yields 



^c-Ar(*'X) 



Am 







ds (^)Q^^(n 



In 



+ 



11 



6K Jo 



2R + sJ 12 
In' ' 



+ 



97 



+ ■ 



16m 

^r' 

1912 m 



315 R 



roo r / 



2R + sJ ' 60 

s 



2R + sJ ^ 6 



(5-3) 



Notice that the 7^-dependent term from the outer integral exactly cancels that from 
the near-zone integral. The other terms, involving integrals over the past history of 
the source are the "tails" , and are in complete agreement with tail terms derived by 
other methods. In addition by suitably combining the tail terms with the leading 
order quadrupole (Q*"'), octopole (Q*-'^), and current quadrupole {J""^) terms in the 
waveform, one can show that part of the effect of these terms, say for a circular 
orbit of frequency uj, is to convert the phase of the wave from ip = oj[t — R) to 
ip = u){t — R — 2mlnR — const). This demonstrates that, despite the use of a flat- 
spacetime wave equation to embody Einstein's equations, the propagation of the 
radiation follows the true, Schwarzschild-like null cones far from the source. 

There are also contributions to the waveform due to intrinsic spin of the bodies, 
which, for compact bodies, occur at 0(e^'^) (spin-orbit) and O(e^) (spin-spin); these 
have been calculated elsewhere. 113) 

Given the gravitational waveform, one can compute the rate at which energy is 
carried off by the radiation (schematically / hhdf2, the gravitational analog of the 
Poynting flux) . For the special case of non-spinning bodies moving on quasi-circular 
orbits (i.e. circular apart from a slow inspiral), the energy flux through 2PN order 
has the form 



~dl 



32 



-V 



m /2927 5 
~ 7 V'336" ^ 1^ 
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3/2 
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/m 
V r 
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293383 
9072 


+ 


380 \" 
9V\ 
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The first term is the quadrupole contribution, the second term is the IPN contri- 
bution, the third term, with the coefficient 47r, is the "tail" contribution, and the 
fourth term is the 2PN contribution first reported jointly by Blanchet et al.c^ 

Similar expressions can be derived for the loss of angular momentum and linear 
momentum. These losses react back on the orbit to circularize it and cause it to 
inspiral. Radiation of linear momentum can also cause "radiation recoil", a phe- 
nomenon which is probably unobservable. The result is that the orbital phase (and 
consequently the gravitational- wave phase) evolves non-linear ly with time. It is the 
sensitivity of the broad-band LIGO and VIRGO-type detectors to phase that makes 
the higher-order contributions to dE/dt so observationally relevant. A ready-to- use 
set of formulae for the 2PN gravitational waveform template, including the non- 
linear evolution of the gravitational- wave frequency (not including spin effects) have 
been published EIP and incorporated into the Gravitational Radiation Analysis and 
Simulation Package (GRASP), a publically available software toolkit. l3) 
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